% File estima_out
% followed after estima_loop
% analyzing reults from estima
% record settings and write them to an excel file. 

% =========================================================================
flag_Hessiandone =0; 
flag_Spectrumdone=0; 
flag_Tablesdone  =0; 
%% 1. Change of path from Xi's computer
cucd = pwd; 
if isunix == 0;
    if ~isempty( findstr(outpath,'home') );
        temp=findstr(upper(outpath),'LJM');
        outpath=cr_dirpath('O:', 'PROJ_LIB', 'LJM',outpath(temp+3:end));
        cd(outpath);
    end
end

% =========================================================================
%% 2. ORDER_DESCEND is the indicator of the ordered modes 
% Creates PARMODEMAT_DESCEND OPTIMAT_DESCEND
cd(cucd); 
logpostdvec=logpostd_mat(:, 2); 
[junk,order_descend]=sort(logpostdvec,'descend'); 
parmodemat_descend=parmodemat(:, order_descend); 
optimat_descend=optimat(:, order_descend); 


save2Priorstrv(parmodemat_descend,Y,sample,Ynames,parnames,funcmod,outpath); 


% =========================================================================
%% 3.1 If ADD2FUNC exists, change function handle to this model after
% checking variances are identical 
if exist('add2func','var')==1 && ~isempty(add2func) 
    dispaj('Switching function handle from ',func2str(funcmod)); 
    dispaj('to  ',[func2str(funcmod),add2func]);
    disp(' ');
    disp('Verifying solutions match...');
    pause(0.25);
    [junk,junk,difvar,difc]=checkmodels_lmj(funcmod,...
        str2func([func2str(funcmod),add2func]),parmodemat_descend(:,1),solveopt,addsol);
    if difvar < 1e-7;
        disp('Solutions identical')
    else
        disp('Solutions do not match')
    end    
    if difc < 1e-10;
        disp('Constants identical');
    else
        disp('Constants do not match');
    end
    clear var1 var2 difvar difc
    funcmod=str2func([func2str(funcmod),add2func]);    
    disp('Changed funcmod handle')        
end 


%===================================================================
%% 4. Compute Hessian 
npoints = 15; 
modname = 'lmj_postmax'; 
parmode = parmodemat_descend(:, 1); 
x = parmodemat_descend(parposest, 1); 
step = 2/npoints; 
try
    addsol.flag_mixfreq = flag_mixfreq;
    addsol.flag_gtv = flag_gtv;
catch
    disp('DEFAULT: Model is not mixed frequency')
    addsol.flag_mixfreq=0;
end
posterioropt.outfold=outpath; 
if exist('Hes','var')==0 
    Hes=[]; 
end 
[junk,Hes,stdcell]=posterior_eval_lmjexog(x,step,npoints,modname, ...
    parmode, Hes,parposest,funcmod,Y,trainvec, prpar, prss, solveopt, addsol, ssposest, parnames, posterioropt);
flag_Hessiandone=1; 
close all;

% ========================================================================
%% 5. Create ALL OUTPUT CELLS 
estima_out_tabs; 

% ========================================================================
%% 6. Ploting STD, Vardecom, IRF, AC, KF state, KF Error for the best mode
parmode_top = parmodemat_descend(:, 1); 
if flag_mixfreq == 0;
    [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames, stanames ,shonames]...
        =feval(funcmod,parmode_top,solveopt,addsol);
else
    [GG, RR, CONS, eu, SDX, ZZ_zero, out, ssvec, flag, ssnames, stanames ,shonames]...
        =feval(funcmod,parmode_top,solveopt,addsol);
    ZZ  =[ZZ_zero;out.AA];
    % added by XI Luo on 03/29/2011
    if ~isfield(out, 'case_tag') % if case_tag does not exist;
        % create one ourselves
        case_tag = nan(out.nzLF, 1);
        if out.flag_tagsum == 0;
            case_tag = 2*ones(out.nzLF, 1); % case 2: average
        elseif out.flag_tagsum == 1;
            for ii = 1:out.nzLF
                if max(out.AA(ii,:)) == 1/out.NUMIT;
                    case_tag(ii,1) = 2;
                elseif max(out.AA(ii,:)) == 1;
                    case_tag(ii,1) = 3;
                end
            end
        end
        out.case_tag = [-1*ones(out.nzHF, 1); case_tag]; % no transformation on high frequency
     end
end
%% 7. Handle Positions of variables to report 
% Check all variables that are to be reported exist
znames=cell(size(ZZ,1),1);
pos_obsinstates=zeros(length(znames),1); 
decompstates_names=cell(size(znames,1),1); 
pos_decompstates=zeros(length(znames),1); 
for ii=1:length(znames)
    temp      =find(ZZ(ii,:)~=0 ) ;
    if length(temp) > 1
        error('Need 1 entry only for each row of ZZ')
    end
    znames(ii)=stanames( temp );
    decompstates_names(ii)=stanames( temp ); 
    pos_obsinstates(ii)=temp; 
    pos_decompstates(ii)=temp; 
    if flag_mixfreq==1 
        if ii > out.nzHF 
            pos_decompstates(ii)=size(ZZ,2)+(ii-out.nzHF); 
            decompstates_names(ii)={ [znames{ii},' TAG'] }; 
        end 
    end 
end
if length(znames)~=size(Y,2); 
    error('Mistmatch ZNAMES and Cols in Y') 
end 
cellposition(state_sel,stanames); 
%cellposition(specrep  ,stanames); 
ana_in.state_sel = state_sel;
ana_in.specrep   = specrep; 
[pos_statesreport,ana_in.flag_adds]=extractlevels(ana_in.state_sel,stanames);

disp('Important Pointers'); 
disp('pos_statesreport=position of states in state_sel'); 
disp('pos_obsinstates =position of observables in states'); 

cd(outpath); 
save workspace; 
cd(cucd); 
disp('Hessian Done'); 


% ========================================================================
%% 8. Spectrum 
ana_in.lowper =6; 
ana_in.highper =32; 
ana_in.npoints =2000; 
figname = ['topmode_specdecom'];
[junk1,spec_outcell,hand_vec]=dsge_spectrum(GG,RR,ZZ,SDX,ana_in.specrep,stanames,shonames,znames,...
    ana_in.lowper,ana_in.highper,ana_in.npoints,save_path, figname);

specdec_mat   =junk1.vfbdecst;
specdec_report=specrep; 
specdec_mat   =specdec_mat( cellposition( specdec_report, ana_in.specrep) , : ); 
cd(outpath);
% save specdec specdec_mat specdec_names;
% cd(cucd);
% display('done!')

cd(cucd); 
clear junk
[~,subfname]=extrsubf(subf,[cucd,'\',root,'\',spec]); 
MakeLaTeX_table(specdec_mat,specdec_report,shonames,'Decomp',['Spectral Decomposition ',subfname(2:end)],['SpDec ',subfname(2:end)],0,outpath);
% elseif flag_mixfreq == 1;
%     if flag_gtv == 1; 
%        ZZ = [ZZ_zero;  out.AA];
%     else
%        ZZ = ZZ_zero;
%     end
%     [outq,outm,spec_outcell]=dsge_spectrum_mixfreq(GG,RR,ZZ,SDX,specrep,stanames,shonames,...
%     znames,ana_in.lowper,ana_in.highper,ana_in.npoints,3,save_path,figname);
% end
spec_outcell=merge_cells(topcell, spec_outcell, 2);
close all
flag_Spectrumdone=1; 
clear junk1 spec_outcell spec_outcell; 



% ========================================================================
%% 8. Theoretical Moments 
cd(cucd)
[vdcom,vcvth,sdevf,ACzero,ACone,vcshocks,auone,vdech,vdech_st,vdcom_st]=vardecom_mode(GG,RR,SDX,ZZ,50);
var_outcell  = crtcell( vdcom, znames , shonames ,'Theoretical Variance Decomposition');
ac0th_outcell= crtcell( ACzero, znames, znames   ,'Theoretical AC(0)'); 
ac1th_outcell= crtcell( ACone , znames, znames   ,'Theoretical AC(1)'); 
cd(cucd); 
stdvsdat_outcell= crtcell([std(Y)' sdevf],znames   ,{'Data Std','Theoretical Std'},'Standard Deviation');  
ac1vsdat_outcell= crtcell([diag(accov(Y,1)) diag(ACone)],znames,{'Data AC(1)','Theoretical AC(1)'},'AC(1)');  

theoretical_cell=merge_cells(topcell,stdvsdat_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,ac1vsdat_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,var_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,ac0th_outcell,2); 
theoretical_cell=merge_cells(theoretical_cell,ac1th_outcell,2); 
clear var_outcell  ac0th_outcell ac1th_outcell stdvsdat_outcell ac1vsdat_outcell; 

par_cell = [parnames, num2cprec(parmode_top, 5)];
ss_cell = [ssnames, num2cprec(ssvec)];
mode_cell = merge_cells(par_cell, ss_cell, 2);
mode_cell = merge_cells(topcell, mode_cell, 2);
mode_cell = merge_cells(mode_cell, stdcell(:, 1:3), 2); 

%return 


%% 9. Save Moments of MODE 
cd(outpath)
if ~isunix
    xlswrite(mode_filename, mode_cell,'Details');
else
    savecell(mode_cell,[],outpath,['Details_', mode_filename]);
end

clear mode_cell; 
try
    if ~isunix
        xlswrite(mode_filename,spec_outcell,'spectrum');
    else
        savecell(spec_outcell,[],outpath,['spectrum','_', mode_filename]);
    end
    disp('Saved Spectrum'); 
    
catch
    disp('Spectrum Tab not saved')
end

try
    if ~isunix 
        xlswrite(mode_filename,theoretical_cell,'TheoreticalMoments');
    else
        savecell(theoretical_cell,[],outpath,['TheoreticalMoments', '_', mode_filename]);
    end
    
    disp('Save Theoretical_Cell'); 
    clear theoretical_cell; 
catch
    disp('Stationary Variance Tab not saved')
end

cd(cucd); 
flag_Tablesdone  =1; 
flag_Complete    =0; 

% cd(outpath); 
% save workspace; 
% cd(cucd); 
close all; 

estima_states; 


cd(outpath); 
save workspace; 
cd(cucd); 
disp('Hessian Done'); 

estima_out_analysis 
